##################################################################
# File:      Replication_Code.R
# Purporse:  Replication of main text figures and tables for 
#                 Public Opinion Quarterly
# Author:    Harry Oppenheimer, Ryan Shandler, Nadiya Kostyuk
# Date:      09/08/22
# Last Edit: 09/12/22
# Data In:   ~/Replication
# Data Out:  ~/Replication/Main Text
# Prev file: None
# Status:    In Progress.
##################################################################


library(psych)
library(tidyverse)
library(cjoint)
library(cregg)

set.seed(02138)

#Set the working directory to the replication files location
setwd("~/Dropbox/Cyber-Terrorism Conjoint Experiment/Replication")
dir.create("Main Text") 


#Functions to make plots from cregg look like the plots from cjoint

make_feature_headers <- function(x, fmt = "%s:") {
  feature_levels <- rev(split(x$level, x$feature))
  for (i in seq_along(feature_levels)) {
    feature_levels[[i]] <- levels(x$level)[match(feature_levels[[i]], levels(x$level))]
    feature_levels[[i]] <- c(feature_levels[[i]], sprintf(fmt, names(feature_levels)[i]))
  }
  factor(as.character(x$level), levels = unique(unname(unlist(feature_levels))))
}

theme_bw1 <- function(base_size = 11, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(axis.text.x = element_text(size = base_size*.9, colour = "black",  hjust = .5 , vjust=1),
          axis.text.y = element_text(size = base_size, colour = "black", hjust = 0 , vjust=.5,family=NULL), 
          axis.ticks = element_line(colour = "grey50"),axis.title.y = element_text(size = base_size,angle=90,vjust=.01,hjust=.1,
                                                                                   family="serif"),plot.title = element_text(face = "bold",family="serif"),legend.position = "none")}

plot.cj_diffs <-
  function(
    x,
    group = attr(x, "by"),
    feature_headers = TRUE,
    header_fmt = "%s:",
    size = 1.0,
    xlab = "Estimated Difference",
    ylab = "",
    legend_title = if (is.null(group)) "Feature" else group,
    legend_pos = "bottom",
    xlim = NULL,
    vline = 0,
    vline_color = "gray",
    theme = theme_bw1(),
    ...
  ) {
    
    # optionally, add gaps between features
    if (isTRUE(feature_headers)) {
      x$level <- make_feature_headers(x, fmt = header_fmt)
      to_merge <- data.frame(feature = unique(x$feature), level = sprintf(header_fmt, unique(x$feature)))
      if ("BY" %in% names(x)) {
        to_merge <- do.call("rbind", lapply(unique(x[["BY"]]), function(lev) {
          to_merge[["BY"]] <- lev
          to_merge
        }))
      } else if (!is.null(group)) {
        to_merge <- do.call("rbind", lapply(unique(x[[group]]), function(lev) {
          to_merge[[group]] <- lev
          to_merge
        }))
      }
      x <- merge(x, to_merge, all = TRUE)
    }
    
    if (is.null(group)) {
      p <- ggplot2::ggplot(data = x, ggplot2::aes_string(x = "estimate", y = "level", colour = "feature"))
    } else {
      p <- ggplot2::ggplot(data = x, ggplot2::aes_string(x = "estimate", y = "level", colour = group, group = group))
    }
    
    if (is.null(xlim)) {
      xmin <- min(x$lower, na.rm = TRUE)
      xmin <- if (xmin < 0) 1.04*xmin else .96*xmin
      xmax <- max(x$upper, na.rm = TRUE)
      xmax <- if (xmax > 0) 1.04*xmax else .96*xmax
      # make symmetric
      if (abs(xmin) > abs(xmax)) {
        xmax <- abs(xmin)
      } else {
        xmin <- -xmax
      }
      xlim <- c(xmin, xmax)
    }
    
    if (!is.null(vline)) {
      p <- p + ggplot2::geom_vline(xintercept = vline, size=.5,colour="black",linetype="dotted")
    }
    
    p <- p + ggplot2::geom_point(position = ggstance::position_dodgev(height = 0.75), size = 3, na.rm = TRUE) +
      ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower", xmax = "upper"),  
                              size = 0.5, height = 0, na.rm = TRUE,
                              position = ggstance::position_dodgev(height = 0.75))
    if (is.null(group)) {
      p <- p + ggplot2::scale_colour_discrete(guide = ggplot2::guide_legend(title = legend_title))
    } else {
      p <- p + ggplot2::scale_colour_discrete(breaks = levels(x[[group]]),
                                              labels = levels(x[[group]]),
                                              guide = ggplot2::guide_legend(title = legend_title))
    }
    p <- p +
      ggplot2::scale_x_continuous(limits = xlim, oob = scales::rescale_none) +
      ggplot2::xlab(xlab) + 
      ggplot2::ylab(ylab) + 
      theme + ggplot2::theme(
        legend.position = "none",
        #  panel.grid.major = ggplot2::element_blank(),
        #  panel.grid.minor = ggplot2::element_blank()
      ) + 
      ggplot2::guides(colour = ggplot2::guide_legend(title = legend_title))
    return(p)
  }


#Function to recode our variables
recode2 <- function(variable, reverse=FALSE, maxVal, minVal, binarize=FALSE, x=NA, x2){
  if (is.factor(variable)){variable <- as.numeric(variable)}
  if (missing(maxVal)){maxVal <- max(variable, na.rm=TRUE)}
  if (missing(minVal)){minVal <- min(variable, na.rm=TRUE)}
  variable[variable > maxVal] <- NA
  variable[variable < minVal] <- NA
  if (reverse){
    temp <- variable
    for (j in 1:maxVal){
      temp[which(variable==j)] <- maxVal-j+1
    }
    return(temp)
  }
  if (binarize){
    if (is.na(x)){stop("Specify the value to dichotomize on")}
    temp <- variable
    if (missing(x2)){
      i <- which(variable == x)
      k <- which(variable < x | variable > x)
    }
    else{
      i <- which(variable >= x & variable <= x2)
      k <- which(variable < x | variable > x2)
    }
    temp[i] <- 1
    temp[k] <- 0
    return(temp)
  }
  return(variable)
}

#Rescale from 0-1
rescale <- function(x){
  return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}

#Convert verbal likert to numeric
likert <- function(x, points=7){
  if (points!=7 & points!=5){stop("Currently requires a 5 or 7 point scale")}
  y <- rep(NA, length(x))
  if (points==7){
    y[which(x=="Strongly agree")] <- 7
    y[which(x=="Agree")] <- 6
    y[which(x=="Slightly agree")] <- 5
    y[which(x=="Somewhat agree")] <- 5
    y[which(x=="Neither agree nor disagree")] <- 4
    y[which(x=="Somewhat disagree")] <- 3
    y[which(x=="Disagree")] <- 2
    y[which(x=="Strongly disagree")] <- 1
  }
  else if (points==5){
    y[which(x=="Strongly agree")] <- 5
    y[which(x=="Agree" | x=="Somewhat agree")] <- 4
    y[which(x=="Neither agree nor disagree")] <- 3
    y[which(x=="Disagree" | x=="Somewhat disagree")] <- 2
    y[which(x=="Strongly disagree")] <- 1	
  }
  return(y)
}

data <- read.csv("Combined Conjoint Data.csv")

### Demographics

data$ID <- 1:nrow(data)

## Gender
data$Male <- ifelse(data$Gender == 2, 1, 0)

## Age
data$age <- 2020 - as.numeric(as.character(data$Year_Born))
data$age_1 <- recode2(data$age, binarize=TRUE, x=18, x2=24)
data$age_2 <- recode2(data$age, binarize=TRUE, x=25, x2=44)
data$age_3 <- recode2(data$age, binarize=TRUE, x=45, x2=64)
data$age_4 <- recode2(data$age, binarize=TRUE, x=65, x2=105)

## Party ID
#data$pid <- NA
#data$pid1 <- rescale(data$Party_ID_US)

#Different coding in the UK and US versions
#data <- data[-which(data$Finished==0),]   # Uncomment for UK version
#data <- data[-which(data$Finished=="False"),]

# ebalance(Treatment, X, base.weight = NULL,
#          norm.constant = NULL, coefs = NULL,
#          max.iterations = 200,
#          constraint.tolerance = 1, print.level = 0)

colnames(data)[87:91] <- c("MethodOfAttack1", "TargetOfAttack1", "ActorType1", "MotivationOfAttack1", "OutcomeOfAttack1")

data_long <- data %>% select(ID, Method1:Outcome7, colnames(data)[grep("Scenario", colnames(data))]) %>%
  gather("Type", "Value", -ID) %>%
  mutate(Round = as.numeric(gsub("[^\\d]+", "", Type, perl = T)), Type = gsub("[0-9]+", "\\1", Type)) %>% filter(grepl("Space", Type)==F) %>%
  spread(Type, Value) %>% mutate(TerrorBinary = ifelse(Scenario == "1", 1, 0)) %>% select(-Scenario)

data_long$Method_short <- dplyr::recode(as.character(data_long$Method), 
                                        "1" = "attack",
                                        "2" = "virus",
                                        "3" = "phishing",
                                        "4" = "USB", 
                                        "5" = "cyberattack")

data_long$Target_short <- dplyr::recode(as.character(data_long$Target), 
                                        "1" = "power station",
                                        "2" = "military facility",
                                        "3" = "shopping mall",      
                                        "4" = "government office building",
                                        "5" = "databases")

data_long$Agent_short <- dplyr::recode(as.character(data_long$Agent), 
                                       "1" = "foreign country",
                                       "2" = "org. w/. foreign government ties",
                                       "3" = "hackers w/. foreign government ties",
                                       "4" = "unaffiliated group of hackers",
                                       "5" = "individual hacker",
                                       "6" = "unknown actor",
                                       "7" = "Jihadi Islamic group")

data_long$Motivation_short <- dplyr::recode(as.character(data_long$Motivation), 
                                            "1" = "unknown",
                                            "2" = "government overthrow",
                                            "3" = "change government policy",
                                            "4" = "revenge")

data_long$Outcome_short <- dplyr::recode(as.character(data_long$Outcome), 
                                         "1" = "no damage",
                                         "2" = "minor explosion",
                                         "3" = "major explosion",
                                         "4" = "theft of millions of dollars",
                                         "5" = "theft and release of\n sensitive information")

firstup <- function(x) {
  substr(x, 1, 1) <- toupper(substr(x, 1, 1))
  x
}

data_long$Method <- firstup(data_long$Method_short)
data_long$Target <- firstup(data_long$Target_short)
data_long$Agent <- firstup(data_long$Agent_short)
data_long$Motivation <- firstup(data_long$Motivation_short)
data_long$Outcome <- firstup(data_long$Outcome_short)

data_long$Method <- as.factor(data_long$Method)
data_long$Target <- as.factor(data_long$Target)
data_long$Agent <- as.factor(data_long$Agent)
data_long$Motivation <- as.factor(data_long$Motivation)
data_long$Outcome <- as.factor(data_long$Outcome)

data_long$Method_short <- as.factor(data_long$Method_short)
data_long$Target_short <- as.factor(data_long$Target_short)
data_long$Agent_short <- as.factor(data_long$Agent_short)
data_long$Motivation_short <- as.factor(data_long$Motivation_short)
data_long$Outcome_short <- as.factor(data_long$Outcome_short)

#attribute_list <- list()
#attribute_list[["Method_short"]] <- as.character(unique(data_long$Method_short))
#attribute_list[["Target_short"]] <- as.character(unique(data_long$Target_short))
#attribute_list[["Agent_short"]] <- as.character(unique(data_long$Agent_short))
#attribute_list[["Motivation_short"]] <- as.character(unique(data_long$Motivation_short))
#attribute_list[["Outcome_short"]] <- as.character(unique(data_long$Outcome_short))

#constraint_list <- list()
## Constraints on Target and Outcome
## Cannot have outcome = "minor explosion" or "major explosion" if target =  "databases"
#constraint_list[[1]] <- list()
#constraint_list[[1]][["Outcome_short"]] <- c("minor explosion", "major explosion")
#constraint_list[[1]][["Target_short"]] <- c("databases")
#
#cyberdesign <- makeDesign(type='constraints', attribute.levels=attribute_list,
#                          constraints=constraint_list)
#

attribute_list <- list()
attribute_list[["Method"]] <- as.character(unique(data_long$Method))
attribute_list[["Target"]] <- as.character(unique(data_long$Target))
attribute_list[["Agent"]] <- as.character(unique(data_long$Agent))
attribute_list[["Motivation"]] <- as.character(unique(data_long$Motivation))
attribute_list[["Outcome"]] <- as.character(unique(data_long$Outcome))

constraint_list <- list()
# Constraints on Target and Outcome
# Cannot have outcome = "minor explosion" or "major explosion" if target =  "databases"
constraint_list[[1]] <- list()
constraint_list[[1]][["Outcome"]] <- c("Minor explosion", "Major explosion")
constraint_list[[1]][["Target"]] <- c("Databases")

cyberdesign <- makeDesign(type='constraints', attribute.levels=attribute_list,
                          constraints=constraint_list)

data$cpu1 <- ifelse(data$CountryCode == 3, 
                    ifelse(data$Comp_Literacy1 == 4|data$Comp_Literacy1 == 2, 1, 0), 
                    ifelse(data$Comp_Literacy1 == 4, 1, 0))

data$cpu2 <- ifelse(data$Comp_literacy2 == 1, 1, 0)
data$cpu3 <- ifelse(data$Comp_literacy3 == 4, 1, 0)

data$cpu1[which(is.na(data$cpu1)==T)] <- 0
data$cpu2[which(is.na(data$cpu2)==T)] <- 0
data$cpu3[which(is.na(data$cpu3)==T)] <- 0

data$cpu_knowledge <- data$cpu1 + data$cpu2 + data$cpu3

data_cj <- data %>% mutate(ideorate_numeric = as.numeric(as.character(Political_position)),
                           liberal = as.numeric((ideorate_numeric==1|ideorate_numeric==2|ideorate_numeric ==3)),
                           conservative = as.numeric((ideorate_numeric==4|ideorate_numeric==5|ideorate_numeric==6)),
                           terror_exp1 = as.numeric(as.character(Prev_Exposure_Terror_1)) -1,
                           terror_exp2 = as.numeric(as.character(Prev_Exposure_Terror_2)) -1,
                           terror_exp_mean = (terror_exp1 + terror_exp2)/2,
                           cyber_exp1  = as.numeric(as.character(Prev_Exposure_Cyber_1)),
                           cyber_exp2  = as.numeric(as.character(Prev_Exposure_Cyber_2)),
                           cyber_exp3  = as.numeric(as.character(Prev_Exposure_Cyber_3)),
                           cyber_exp_mean = (cyber_exp1 + cyber_exp2 + cyber_exp3)/3) %>%
  select(ID, Male, liberal, conservative, terror_exp_mean, cyber_exp_mean, Education, CountryCode) %>%
  right_join(data_long, by = "ID")

data_cj$CountryCode <- factor(data_cj$CountryCode)

data_cj$Motivation <- factor(data_cj$Motivation, levels = c("Revenge", "Government overthrow", "Change government policy", "Unknown"))
data_cj$Outcome <-   factor(data_cj$Outcome, levels = c("No damage", "Major explosion", "Theft and release of\n sensitive information",
                                                        "Minor explosion", "Theft of millions of dollars"))
data_cj$Method <-    factor(data_cj$Method, levels = c("Cyberattack", "Attack", "USB", "Phishing", "Virus"))
data_cj$Target <-    factor(data_cj$Target, levels = c("Military facility", "Databases", "Government office building", "Power station", "Shopping mall"))
data_cj$Agent <-     factor(data_cj$Agent, levels =  c("Foreign country","Hackers w/. foreign government ties", "Jihadi Islamic group",
                                             "Org. w/. foreign government ties","Unaffiliated group of hackers",
                                             "Individual hacker", "Unknown actor"))

#############
## Paper Figure 2
#############

results_pooled1 <- cjoint::amce(TerrorBinary ~ Method + Target + Agent +  Motivation + Outcome, data=data_cj,
                                cluster=TRUE, respondent.id="ID", design=cyberdesign)

plot(results_pooled1, xlim = c(-.15, .15), xlab = "Change in Pr(cyberterrorism)")

ggsave("Main Text/figure2.tiff", width = 10, height = 8)


#############
## Paper Figure 3
#############
data_cj$Country <- NA
data_cj$Country[which(data_cj$CountryCode == 1)] <- "United States"
data_cj$Country[which(data_cj$CountryCode == 2)] <- "United Kingdom"
data_cj$Country[which(data_cj$CountryCode == 3)] <- "Israel"
data_cj$Country <- factor(data_cj$Country, levels = c("Israel", "United States", "United Kingdom"))

results_pooled_country <- cjoint::amce(TerrorBinary ~ Method:Country + Target:Country + Agent:Country +
                                 Motivation:Country + Outcome:Country, data=data_cj,
                               cluster=TRUE, respondent.id="ID",  respondent.varying = "Country", design=cyberdesign)

plot(results_pooled_country, xlim = c(-.2, .2),
     xlab = "Change in Pr(cyberterrorism)",plot.display="interaction", 
     group.order=c("Agent","Method","Motivation",
                   "Outcome","Target"))

ggsave("Main Text/figure3.tiff", width = 10, height = 8)


#############
##Paper Figure 4
#############

data_cj$CountryCode <- factor(data_cj$CountryCode)
data_cj$CountryName <- recode(data_cj$CountryCode, "1" = "US", "2" = "UK", "3" = "IL")
data_cj$CountryName_2 <- factor(data_cj$CountryName, levels = c("UK", "IL", "US"))


results_interact_diff <- cregg::cj(TerrorBinary ~  Method  + Agent +  Motivation + Target * Outcome,
                            data=data_cj, cluster=TRUE, respondent.id="ID", by = ~CountryName,
                            feature_order = c("Agent", "Method", "Motivation", "Outcome", "Target"),
                            estimate = "mm_diff")

results_interact_diff$feature <- as.character(results_interact_diff$feature)
results_interact_diff$feature <- gsub("\\_.*","",results_interact_diff$feature)

ukus <- results_interact_diff[grep("UK - US", results_interact_diff$BY),]
ilus <- results_interact_diff[grep("IL - US", results_interact_diff$BY),]

ukus <- ukus[c(which(ukus$level == "Unknown actor"),
               which(ukus$level == "Individual hacker"),
               which(ukus$level == "Unaffiliated group of hackers") ,
               which(ukus$level == "Org. w/. foreign government ties"),
               which(ukus$level == "Jihadi Islamic group"),
               which(ukus$level == "Hackers w/. foreign government ties"),
               which(ukus$level == "Foreign country"),
               which(ukus$level == "Virus"),
               which(ukus$level == "Phishing"),
               which(ukus$level == "USB"),
               which(ukus$level == "Attack"),
               which(ukus$level == "Cyberattack"),
               which(ukus$level == "Unknown"),
               which(ukus$level == "Change government policy"),
               which(ukus$level == "Government overthrow"),
               which(ukus$level == "Revenge"),
               which(ukus$level == "Theft of millions of dollars"),
               which(ukus$level == "Minor explosion"),
               which(ukus$level == "Theft and release of\n sensitive information"),
               which(ukus$level == "Major explosion"),
               which(ukus$level == "No damage"),
               which(ukus$level == "Shopping mall"),
               which(ukus$level == "Power station"),
               which(ukus$level == "Government office building"),
               which(ukus$level == "Databases"),
               which(ukus$level == "Military facility")) ,]

ilus <- ilus[c(which(ilus$level == "Unknown actor"),
               which(ilus$level == "Individual hacker"),
               which(ilus$level == "Unaffiliated group of hackers") ,
               which(ilus$level == "Org. w/. foreign government ties"),
               which(ilus$level == "Jihadi Islamic group"),
               which(ilus$level == "Hackers w/. foreign government ties"),
               which(ilus$level == "Foreign country"),
               which(ilus$level == "Virus"),
               which(ilus$level == "Phishing"),
               which(ilus$level == "USB"),
               which(ilus$level == "Attack"),
               which(ilus$level == "Cyberattack"),
               which(ilus$level == "Unknown"),
               which(ilus$level == "Change government policy"),
               which(ilus$level == "Government overthrow"),
               which(ilus$level == "Revenge"),
               which(ilus$level == "Theft of millions of dollars"),
               which(ilus$level == "Minor explosion"),
               which(ilus$level == "Theft and release of\n sensitive information"),
               which(ilus$level == "Major explosion"),
               which(ilus$level == "No damage"),
               which(ilus$level == "Shopping mall"),
               which(ilus$level == "Power station"),
               which(ilus$level == "Government office building"),
               which(ilus$level == "Databases"),
               which(ilus$level == "Military facility")) ,]

results_interact_diff_2 <- cregg::cj(TerrorBinary ~  Method  + Agent +  Motivation + Target * Outcome,
                              data=data_cj, cluster=TRUE, respondent.id="ID", by = ~CountryName_2,
                              feature_order = c("Agent", "Method", "Motivation", "Outcome", "Target"),
                              estimate = "mm_diff")
results_interact_diff_2 <- results_interact_diff_2[grep("IL - UK", results_interact_diff_2$BY),]

results_interact_diff_2$feature <- as.character(results_interact_diff_2$feature)
results_interact_diff_2$feature <- gsub("\\_.*","",results_interact_diff_2$feature)

results_interact_diff_2 <- results_interact_diff_2[c(which(results_interact_diff_2$level == "Unknown actor"),
               which(results_interact_diff_2$level == "Individual hacker"),
               which(results_interact_diff_2$level == "Unaffiliated group of hackers") ,
               which(results_interact_diff_2$level == "Org. w/. foreign government ties"),
               which(results_interact_diff_2$level == "Jihadi Islamic group"),
               which(results_interact_diff_2$level == "Hackers w/. foreign government ties"),
               which(results_interact_diff_2$level == "Foreign country"),
               which(results_interact_diff_2$level == "Virus"),
               which(results_interact_diff_2$level == "Phishing"),
               which(results_interact_diff_2$level == "USB"),
               which(results_interact_diff_2$level == "Attack"),
               which(results_interact_diff_2$level == "Cyberattack"),
               which(results_interact_diff_2$level == "Unknown"),
               which(results_interact_diff_2$level == "Change government policy"),
               which(results_interact_diff_2$level == "Government overthrow"),
               which(results_interact_diff_2$level == "Revenge"),
               which(results_interact_diff_2$level == "Theft of millions of dollars"),
               which(results_interact_diff_2$level == "Minor explosion"),
               which(results_interact_diff_2$level == "Theft and release of\n sensitive information"),
               which(results_interact_diff_2$level == "Major explosion"),
               which(results_interact_diff_2$level == "No damage"),
               which(results_interact_diff_2$level == "Shopping mall"),
               which(results_interact_diff_2$level == "Power station"),
               which(results_interact_diff_2$level == "Government office building"),
               which(results_interact_diff_2$level == "Databases"),
               which(results_interact_diff_2$level == "Military facility")) ,]
colnames(results_interact_diff_2)[12] <- "CountryName"


ukus$level <- factor(paste(" ", ukus$level), levels = paste(" ", ukus$level))
ilus$level <- factor(paste(" ", ilus$level), levels = paste(" ", ilus$level))
results_interact_diff_2$level <- factor(paste(" ", results_interact_diff_2$level), levels = paste(" ", results_interact_diff_2$level))

plot.cj_diffs(rbind(ukus, ilus,results_interact_diff_2), 
              vline = 0, xlim = c(-.15, .15)) + ggplot2::facet_wrap(~BY) 

ggsave("Main Text/figure4.tiff", width = 10, height = 8)

